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Abstract 

We examine a phase transition in a model of random spatial permutations which originates in a study 
of the interacting Bose gas. Permutations are weighted according to point positions; the low-temperature 
onset of the appearance of arbitrarily long cycles is connected to the phase transition of Bose-Einstein 
condensates. In our simplified model, point positions are held fixed on the fully occupied cubic lattice 
and interactions are expressed as Ewens-type weights on cycle lengths of permutations. The critical 
temperature of the transition to long cycles depends on an interaction-strength parameter a. For weak 
interactions, the shift in critical temperature is expected to be linear in a with constant of linearity c. 
Using Markov chain Monte Carlo methods and finite-size scaling, we find c = 0.618 ±0.086. This finding 
matches a similar analytical result of Ueltschi and Betz. We also examine the mean longest cycle length 
as a fraction of the number of sites in long cycles, recovering an earlier result of Shepp and Lloyd for 
non-spatial permutations. 



1 Introduction 

The model of random spatial permutations arises in the study of the Bose gas. Its history includes Bose- 
Einstein, Feynman |Feynman] , Penrose- Onsager [POj . Siito |Siitoll ISiito2) . and Betz-Ueltschi-Gandolfo- 
Ruiz [GRU1 IU071 IB UP 71 IBU08] . Such random permutations arise physically when one symmetrizes the 
iV-boson Hamiltonian with pair interactions, then applies a multi-particle Feynman-Kac formula and a 
cluster expansion BU07, BU08 . Specifically, given points xi, . . . , x^ in the box [0, L] 3 and temperature T, 
permutations 7r are given probability weights proportional to the Gibbs factor e~ H ^ where 

T n 

F W = jEll^-^wlli (1.1) 

i=l 

for the non- interacting case. The notation || • ||a indicates the natural distance on the 3-torus: 

||x-y|| A = min{||x-y + in||} (1.2) 

(The sum in equation Ijl.ljl is scaled by temperature rather than reciprocal temperature. This surprising 
feature is opposite that of many models in statistical mechanics.) These energy terms involve lengths of 
permutation jumps; additional interaction terms take the form 

). (1-3) 
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i.e. permutation jumps from sites Xj and Xj interact pairwise. In the above-cited papers of Betz and Ueltschi, 
these may be approximated and rearranged such that one obtains interaction terms of the form 

N 

J2a e r e {n) (1.4) 

1=1 

where 7^(71") counts the number of ^-cycles of the permutation tt, and the coefficients cti are cycle weights. 

When the Bose gas is cooled below a critical temperature T c , there is a phase transition: a macroscopic 
fraction of the bosons are found in the ground state of the external potential, and such particles are said to 
participate in a Bose-Einstein condensate. In the permutation representation, this transition manifests itself 
as the onset of long permutation cycles. Bose suggested the statistics carrying his name for describing the gas 
of photons; Einstein developed the notion of what we now call Bose-Einstein condensation, and computed 
the critical temperature for the non-interacting Bose gas. The critical temperature for liquid helium, where 
intcrparticle interactions are strong, is lower than would be expected BBHLV for non-interacting atoms 
of the same density. For weakly interacting systems, however, an emerging consensus is that interactions 
increase the critical temperature. See BBHLV] [SU09] for surveys. Concretely, for interactions parameterized 
by some a, one defines 

AT c (a) = y ' 



It is well accepted that 



r c (o) 



lim ATJa) = cp 1/3 a 

a— >0 



where p is the particle density, i.e. that for small a the shift in critical temperature is linear in a. What is 
more contentious, as enumerated in the surveys cited above, is the value of the constant c. 

The interaction terms (equation (11.31) ) for the permutation representation of the Bose gas are difficult to 
compute. Moreover, it is interesting to consider the model of random spatial permutations for its own sake. 
In [GRU] , a simulational approach is taken for points held fixed on the fully occupied unit lattice in the non- 
interacting case. In the papers |U071 IBU07] . Betz and Ueltschi examine the Bose-gas permutation weights 
with point positions allowed to vary on the continuum; an exact expression for the critical temperature is 
stated and proved for a simplified interaction model in which only two-cycles interact. That is, interactions 
are of the form of 11.41 with ct2 — a, where a is related to a hard-core scattering length, and the remaining 
cycle weights are zero. In [BU08] . this approach is extended to a model in which all the ae's may vary, but 
with the hypothesis that on goes to zero faster than 1/ log(€). In this paper, we take a simulational approach 
to points on the fully occupied unit lattice, with cycle weights constant in I — removing the decaying-cycle- 
weight hypothesis. The shift in critical temperature is nonetheless found to match that predicted by Betz 
and Ueltschi. 

An outline of the paper is as follows. Section [5] provides background necessary to understand the results 
of the paper: the probability model is defined in section 12.11 qualitative and quantitative behavior of long 
cycles are discussed in sections 12.21 and 12.31 respectively. Known results and conjectures are listed in section 
12.41 In section [3] the simulational methods are presented. The swap-only and swap-and-reverse algorithms 
generate simulational data; these algorithms are proved correct in sections 13.11 through 13.31 The finite-size- 
scaling method, which reduces the raw simulational data, is summarized in section T3. 5 1 Section SI presents 
the data and its analysis in full detail: estimation of critical exponents and critical temperature in sections 
14.11 through 14.31 verification of the finite-size-scaling hypothesis in section 14.41 and final results in sections 
S3] through SJl 
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Figure 1: A spatial permutation. 



2 The model of random spatial permutations 

Here we review concepts from [BU07, BU08 , fixing notation and intuition to be used in the rest of the paper. 
2.1 The probability model 

The state space is Qa,n — A N x Sn, where A = [0, L} 3 with periodic boundary conditions; point positions 
are X = (xi, . . . , xat) for Xi, . . . , £ A. The Hamiltonian takes one of two forms. In the first, relevant to 
the Bose gas, we have 



where T = 1//3 and the V terms are interactions between permutation jumps. (The temperature scale factor 
T/4, not /3/4, is surprising but correct for the Bose-gas derivation of the Hamiltonian.) In the second form 
of the Hamiltonian, considered in this paper, we use interactions which are dependent solely on cycle length: 



where re(Tt) is the number of ^-cycles in it and the a/s are free parameters, called cycle weights. One 
ultimately hopes to choose the a/s appropriately for the Bose gas; even if not, the model is well-defined and 
of its own interest. 

Different choices of ai result in different models: The non-interacting model |GRUj has at = 0. The two- 
cycle model [BU071 [XJ07] . has a2 = a and other cycle weights are zero. The general-cycle model has no 
restrictions on ai. In [BU08] , the decaying cycle-weight case of the general-cycle model is considered: the 
only restriction on at is that ai goes to zero in t faster than 1 / log £. The Ewens model, treated in this paper 
(see also [Ewens ), is another special case of the general-cycle model: it has ag = a constant in I. 

One may hold point positions fixed, e.g. on the fully occupied unit lattice; this approach has been taken for 
all simulations done up to the present by the author and by Gandolfo [GRU , including specifically the work 
described in this paper. One obtains a Gibbs probability distribution on Sn- 




(2.1) 




(2.2) 




e 



H(X,7r) 



F(A,X) ' 



(2.3) 
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(Alternatively, one may integrate over all positions in A, with a resulting Gibbs distribution on Sn. Here, 
several analytical results are available |BU071 [BU08] .) For a random variable S(tt), we have 

E[S] = J2 PWSfr). (2-4) 
2.2 Qualitative characterization of long cycles 

One next inquires which permutations are typical in this temperature-dependent probability distribution on 
Sn. In this section we develop intuition; in the next section, we construct quantitative descriptions of the 
ideas presented here. 




Figure 2: Some typical permutations for high T, medium but subcritical T, and low T. 

As T —¥ oo, the probability measure becomes supported only on the identity permutation: the distance- 
dependent terms are large whenever any jump has non-zero length. For large but finite T, there are tiny 
islands of 2-cycles, 3-cycles, etc. On the other hand, as T — > 0, length-dependent terms go to zero, and 
the probability measure approaches the uniform distribution on Sn- the distance-dependent terms all go 
to zero. For intermediate T, one observes that the length ||7r(x) — x||a of each permutation jump remains 
small, increasing smoothly as T drops. 

For T above a critical temperature T c , all cycles are short: two-cycles, three-cycles, and so on. We find 
T c 6.86 at a = 0, and positive a terms increase T c . At T c , though, there is a phase transition: for T <T C 
jump lengths remain short but long cycles form. Quantitatively, let l max be the length of the longest cycle in 
7r, with E[^ max ] its mean over all permutations. We observe that for T > T c , E[£ max ] grows only perhaps as 
fast as log(i) as L — > oo. That is to say, E[£ max ]/N goes to zero as N — > oo. For T < T c , on the other hand, 
E[4nax] scales with N, i.e. E[£ max ]/A approaches a temperature-dependent constant as N — > oo: there are 
arbitrarily long cycles, or infinite cycles, in the infinite- volume limit. See figure [2] for depictions of typical 
permutations at high T, subcritical T, and low T; see figure [3] for plots of K[£ max ]/N as a function of T for 
various system sizes with N — L 3 . Note in particular that higher alpha shifts the order-parameter curve to 
the right, with resulting upward shift in critical temperature T c . 

Feynman's claim for the Bose gas is that Bose-Einstein condensation occurs if and only if there are infinite 
cycles in the infinite- volume limit. The central point of this approach is that the system energy has been recast 
in terms of permutations, which are amenable to analysis and simulation. This permits a new perspective 
on the venerable question: how does the critical temperature of Bose-Einstein condensation depend on 
inter-particle interaction strength? 

Obtaining a full answer to this notoriously difficult question is a long-term project. As an intermediate step, 
we here consider the Ewens cycle-weight Hamiltonian with point positions on the unit fully occupied unit 
lattice. Through careful use of MCMC algorithms, statistical analysis, and finite-size scaling, we are able to 
quantify the dependence of critical temperature on interaction strength. 
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Figure 3: Order parameter fu = E[£ max ]/iV for finite systems, with a = 0,0.004. Interactions increase the 
critical temperature. The shift is slight, but visible; we work in the regime of small interaction parameters. 
See section H] for a quanitative analysis of this shift. 



2.3 Quantitative characterization of long cycles 

Various order parameters may be defined; all of them may be used to locate the critical temperature T c (a). 
The fraction E[^ max ]/./V discussed above will, for brevity, be hereafter referred to as fu- The fraction of sites 
in long cycles, //, is described in detail in |GRU] . The correlation length £(T) is defined to be the spatial 
length of the cycle containing a given point x: for T < T c , it blows up as L increases. Namely, we define 

1 N 

Sx(tt) = ||tt(x) -x||a and s(tt) = — ^ s Xi (tt). 

i=l 

The expectation over all tt of s x is the same as s, of course; in a Monte Carlo simulation, however, the latter 
yields a larger sample size and thus a smaller error bar. We use £(T) = E[s]. 

Winding numbers count the integer number of x, y, z wraps around the 3-torus (A with periodic boundary 
conditions). Specifically, the winding number of a permutation n is the triple 

1 N 

W = (W x ,W y ,W z ) = -^d A (x, (l) ,x 2 ), (2.5) 

i=i 

where d\ is the difference vector defined as follows. For z € A, we define a zero-centered modulus vector 
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m.i(z). For x, y € A, this gives rise to a difference vector cLa(x, y): 



We also write 



m L {zi) 

m L (z) = ( m L {z 2 ) j (2.6) 

iii(z) = n € Z which minimizes |z + nl/| (2-7) 

m L (z) = z + n L (z)L (2.8) 

d A (x,y) =m A (x-y). (2.9) 

7"2 "117- -117- TJ/2 | tt/2 | ti/2 



w 2 = w • w = wi + Wy + Wf 



The scaled winding number |PC87] is 



, (W 2 )L 2 

is = 



3/3N 

Lastly, the order parameter fw is the fraction of sites which participate in winding cycles. 

The order parameters fi(T), fs{T), and fw(T) show behavior similar to ju '■= E[£ max ]/iV (figure [3]): 
asymptotically as N — > oo, they are zero for T > T c and non-zero for T < T c . For finite N, the curves 
remain analytic: finite-size effects persist. The inverse correlation length 1/£(T), on the other hand, is zero 
for T <T C and non-zero for T > T c . 

Our goal is to quantify the dependence of T c on a, where 

AT c (a) = T ^)~ T ^) . (2 .i0) 
Known results and conjectures are formulated quantitatively in terms of lim Q ^ AT c (a). 



2.4 Known results and conjectures 

Known results for point locations averaged over the continuum are obtained largely using Fourier methods 
BU08 , which are unavailable for point positions held fixed on the lattice. Betz and Ueltschi have determined 
AT c (a), to first order in a, for two-cycle interactions |BU07| and decaying cycle weights [BU08| . (This taps 
into a long and controversial history in the physics literature: see [BBHLV] or [SU09] for surveys.) The 
critical (p,T,a) manifold relates p c to T c . Specifically, 

Pc ( ai ,a 2 , ...) = £ / e-^fW dk = —L_ ^ e -«w^3/ a (2 . U) 

AT c (a) = cp^a, for a w 0. (2.12) 
Using this formula for constant cycle weights at = a and for lattice density p = 1, we have 



(W)3/2 ' ic " C(3/2) 2 / 3 ' 

= ^) 7 r c (0) = e2Q/3 _ 2a 
v ; T c (0) 3 ' 



(2.13) 



We inquire whether this result, obtained for decaying cycle weights with point positions varying on the 
continuum, holds for Ewens weights with point positions held fixed on the lattice. 



G 



For on = (the non-interacting model), M[t max ]/Nfi is constant for T below but near T c . (That is, the two 
order parameters fj and E[^ max ]/iV have the same critical exponent.) For uniform-random permutations 
(Shepp and Lloyd 1966 [EG] solved Golomb's 1964 question jGolombj ). E[£ max ]/iV w 0.6243; unpublished 
work of Betz and Ueltschi has found E[^ max ]/iV// is that same number for the non-interacting case on = 0. 
The intuition is that long cycles are uniformly distributed within the zero Fourier mode. (This was proved 
in section 5 of |Sutol| . Other results on the distribution of the length and number of cycles for probabilities 
depending only on the conjugacy class can be found in sections 2 and 5 of [Siitolj . and in |Siit62j .) We 
conjecture that E[£ max ]/iV// is a-dependent but constant in T (for T below but near T c ) for all interaction 
models. 

We suspect that the fine details of point positions are unimportant for the shift in critical temperature. 
Thus, AT c (a) on the lattice should be similar to that on the continuum, if decaying cycle weights are used. 
For Ewens interactions, though, AT c (a) is theoretically unknown for Ewens interactions with points either 
on the continuum or on the lattice. The simulational treatment in this paper is the only known attack on 
this question. 



3 Simulational methods 

We run Markov chain Monte Carlo experiments for various values of L, T, and interaction strength a. For 
each parameter combination, we generate M typical permutations tt\ , . . . , ttm from the stationary distribu- 
tion, using MCMC algorithms described below, and we compute random variables Xi = X{-Ki). (The values 
of M used are 10 5 away from T c , and 10 6 near T c where sample variance is higher.) We find the sample mean 
and estimate the variance of the sample mean. The correlation of the X^s complicates the latter. Finite-size 
scaling compensates for finite-size effects: mathematically, we are interested in estimating infinite-volume 
quantities based on finite- volume numerical experiments. 

3.1 The swap-only algorithm 

Recall from section |2~T1 that the expectation of a random variable S (such as £, /m, fw, fi, fs) is 

E[S] = ]T P(n)S(w). 

The number of permutations, Nl, grows intractably in N. As is typical in Markov chain Monte Carlo methods 
|Berg| ILB] . one contents oneself with a smaller number of samples: the expectation is instead estimated by 
summing over some number M (10 5 or 10 6 ) of typical permutations. 



x y x y 

• • swap w.p. •> • 

| i ~lAe-- ^ 

tt( x ) f(y) tt(x) ?r(y) 



Figure 4: Metropolis moves for the swap-only algorithm. 

The swap-only algorithm for sampling from the Gibbs distribution (equation (|2.3[) ) is as follows: 

• Start with the identity or uniform-random permutation. 

• Sweep through sites x of the lattice in either lexical or uniform-random order. 



• For each site x, do a Metropolis step: 

— Choose a site 7r(y) from among the six nearest neighbors of 7r(x). 

— Propose to change tt to the permutation tt' which has tt'(z) = ir(z) for all z ^x,y but 7r'(x) = 7r(y) 
and 7r'(y) = 7r(x). (See figure SJ) 

— With probability proportional to min {l,e' AH } where AH = H(tt') - H(ir), accept the change. 
(If the change is rejected, tt' — tt.) 

• After each sweep, obtain a value of each random variable for inclusion in computation of its sample 
mean. 

One starts accumulating data only after a suitable number of thermalization sweeps. The idea is that the 
initial, identity permutation is not typical, nor are the first few afterward. The integrated autocorrelation 
time |Berg| of system energy H gives an idea of how many Metropolis sweeps should be discarded before the 
permutations become typical. Also, one may examine H to ensure that it has reached its long-term average 
value. This is explained in detail in [Kerl]. We next prove correctness of this algorithm. 

3.2 Explicit construction of the Markov matrix 

For section T3. 31 we will need an explicit construction of the Markov matrix corresponding to the swap-only 
algorithm as described in section 13.11 The Markov perspective on the algorithm is that the distribution 
P^(tt) of the first permutation is either supported solely on the identity, or uniform on all A! permutations. 
The distribution for subsequent permutations is 

p( fc+1 V)= P {k) (K)M{n^') 

or, in matrix/vector notation, 

p(fe+i) _ pW]y[_ 

In this section we precisely describe the matrix M; in section |3~31 we show that approaches the Gibbs 
distribution (equation (|2.3[) h 

The matrix M is A! x AH: rows are indexed by m, . . . , 7Tjv! and columns are indexed by tt[, . . . , ir' N ,. Most 
of the entries of M are zero: Metropolis steps change only two permutation sites whereas most tt, tt' differ 
at more than two sites. 

Definition 3.1. For tt,tt' £ Sn, define 

d( 7 r, 7 r') = #{i = l,2,...,iV:7r(i)^ 7 r'(i)}. 

Remark. Note that d(TT,Tr') 7^ 1 since if two permutations agree on A — 1 sites, they must agree on the 
remaining site. It is easily shown that the function (L(tt, it') is a metric on Sn- 

Definition 3.2. Lattice sites x, y are nearest-neighbor if ||x — y|| a = 1- 

Definition 3.3. For tt £ Sn, define 

R(tt) = {tt' e S N : o](tt, tt') = 2 and ||vr(x) - 7r(y)|| A = 1} 

where the x and y are taken to be the two points at which tt, tt' differ. Then R{tt) is the set of permutations 
tt' reachable from tt on a swap. 
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We construct the Markov matrix for use when sites x are selected at uniform random. (The matrices for use 
when x is selected sequentially are similar.) For each tt G Sn, 



-Mia e - H ( 7r ')+^W 

37V 1 



M{tt,tt') = < 



tt' G R{tt), 



(3.4) 



7r"efl(Tr) 



0, 



otherwise. 



To justify the choice of prefactor 1/3N, note that there are N choices of lattice points x. For each x, there 
are 6 choices of 7r(y) which are nearest neighbors to 7r(x). This double-counts the 3N distinct choices of tt' 
reachable from tt in a single Metropolis step, since choosing x and then y results in the same Metropolis 
step as choosing y and then x. 



3.3 Correctness of the swap-only algorithm 

It is clear that the swap-only algorithm produces a sequence of permutations, but with what distribution? 
From Markov-chain theory, we know the following: If the chain is irreducible, aperiodic, and satisfies detailed 
balance, then the chain has the Gibbs distribution (equation (|2.3p ) as its unique invariant distribution. 

We note the following terminology: detailed balance is the same as reversibility . Also, an irreducible, aperiodic 
chain on a finite state space is called ergodic. Also note from Markov-chain theory that all states in a 
recurrence class have the same period. Thus, if we can show that the chain is irreducible (i.e. the entire 
state space is a single recurrence class), then for aperiodicity of the chain it suffices to show that a single 
state (e.g. the identity permutation) has period 1. 

Proposition 3.5 (Irreducibility) . For all tt,tt', there is an n such that M"(7r,7r') > 0. That is, any 
permutation is reachable from any other. 

Proof. Transpositions generate 5/v : for all tt G 5/v, there exist transpositions o~i, . . . ,o~ m such that tt = 
Y[j—i Thus, it suffices to show that given any permutation tt and any two points x and z, so tt : x t— > 7r(x) 
and tt : z H> tt(z), we can construct a sequence of swaps sending tt to tt' so that tt' : x H> tt(z), tt' : z H> 7r(x), 
and 7r'(y) = 7r(y) for all y ^ x, z. (If 7r(x) and 7r(z) are nearest-neighbor lattice sites, of course, then a single 
swap does the job.) 

Define G a .b : Sn —> Sn to be the swap operator for nearest-neighbor lattice sites 7r(a) and 7r(b). Write 
tt' = G^tt. Given x and z, there is a (non-unique) sequence of lattice sites yo,yi,y2> • • • >y™ such that 
yo = x, y„ = z, and ||7r(yi+i) — v(yi)|| A = 1 for i = 0, 1, . . . , n — 1. (See figure [3]) We will construct 
a sequence of swaps along this nearest-neighbor path whose end result is to swap the permutation arrows 
starting at x and z, leaving all other arrows unchanged. We first need a lemma about compositions of swaps. 

Notation 3.6. Given Xi, . . . , x^r and a permutation 7r, we may write tt as an image map with the x^'s along 
the top row and their images along the bottom row: 

( Xi ... x w \ 

y 7r(xi) . . . 7r(xAr) J 

We find that the composition of maps 

(Gy n ,yi ° G y„,y 2 ° • ■ • ° Gy„, y „_ 2 ° Gy„,y„-i) ° (G yo , y „ ° G yn , y „_ 1 o . . . o G yo , y2 o G yo , yi ) (3.7) 
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y — 7r(y ) (1) 
yi — 7r(yi) 
y 2 — 7r(y 2 ) 
y 3 — ?r(y 3 ) 


yo -w- ""(yo) (2) 

yi • / *« T(yi) 

y 2 — ?r(y 2 ) 
y 3 — 7r(y 3 ) 


yo w- ""(yo) (3) 

yi Ay- 7r(yi) 
y 2 ./ V 7r(y 2 ) 
y 3 — „ 7r(y 3 ) 


yo 7r(yo) (4) 

yi v ^(yi) 

y2 • / \« 7r(y 2 ) 
y 3 ./ \ 7r(y 3 ) 


yo \j' 7r(yo) (5) 

yi -K f - T(yi) 
y2 -X> 7r(y2) 

y 3 ./ \ 7r(y 3 ) 


yo •> r ?r(yo) (6) 

yi ""(yi) 

y 2 .-A. 7r(y 2 ) 

y3 • 7r(y3) 



Figure 5: A sequence of (nearest-neighbor) swaps which results in a non- nearest- neighbor swap. 



swaps the images of x = yo and z = y„ while leaving all other images unchanged, that is, 

( yo yi • ■ • y«-i y« \ / yo yi • ■ • y«-i y« \ 
V ^(yo) ""(yi) •■• T(yn-i) ""(yn) / \ »(y») ^(yi) •■• ^(yn-i) i"(yo) )' 

□ 

Remark. Below we will discuss winding numbers, and the empirical fact that the swap-only algorithm 
changes them only rarely. The chain is irreducible but various non-zero transition probabilities can still be 
very small. 

Definition 3.8. The period of tt is 

P(tt) = gcd{n : P(U n = tt | n = u) > 0} 

where Hk is the random variable which is the permutation appearing at the fcth step of the Markov chain. 
We say that tt has period p if it reappears with probability 1 after every p steps. A permutation tt is aperiodic 
if p{tt) = 1. The chain is aperiodic if p(ir) — 1 for every tt. 

Proposition 3.9 (Aperiodicity). The swap-only algorithm's Markov chain is aperiodic. 

Proof. This follows from irreducibility, which says in particular that for every tt, there is an integer m such 
that M m (7r, tt) > 0. Then M n (n, tt) > for all n > m, implying p(ir) = 1. □ 

Proposition 3.10 (Detailed balance). For all tt,tt' £ Sn, 

P(tt)M(tt, tt') = P(tt')M(tt', tt). (3.11) 

Proof. The detailed-balance statement in terms of the Gibbs distribution (equation (|2.3[1 ) and the Metropolis 
transition matrix (equation (|3.4l0 is 

(l A e- H ^e H ^) = e —^- (l A e -*« e W) . 

The Z's cancel. The lemma below shows that M(tt, tt') ^ iff M(tt', tt) ^ 0. If M(tt, tt') = 0, then detailed 
balance holds. If M(n, tt') ^ 0, then there are two cases. If H(tt') < H(tt), then 

e- H ^ (1) = e~ H ^ (e-^)^) . 



10 



If H(n') > H(tt), 



e -ffW (e- H ('') e H W) = e- ff ("') (1). 

In all cases, detailed balance holds. □ 
Lemma 3.12. For all tt, tt' £ Sn, 

M(tt,tt')^0 ^f=> M(tt',tt)^0. 

Proof. This is true since tt' £ R(k) if and only if tt £ R(tt'), which is a direct consequence of the definition 
Elof i?(7r). □ 

This lemma completes the proof that the swap-only algorithm satisfies detailed balance and thus has the 
Gibbs distribution as its invariant distribution. It is not hard to show that if swaps sites x ^ y are in the 
same cycle before a swap, they are in different cycles after the swap, and vice versa. This is not a correctness 
result, but rather a sanity check: it shows that cycles may grow or shrink upon swap-only moves. 




Figure 6: Swaps merge disjoint cycles and split single cycles. The left-hand permutation can be reached 
from the right-hand permutation via a swap, and vice versa. 



3.4 Winding cycles and the swap- and- reverse algorithm 

The propositions of section T3.3I showed that the swap-only algorithm is correct — in particular, any permu- 
tation is reachable from any other with non-zero probability. However, in practice some of these non-zero 
transition probabilities can be quite small. In particular, we observe that the swap-only algorithm almost 
always generates permutations with zero winding number. 

This problem, and a partial solution, is explained intuitively by figure [7] and rigorously in |Kerl] . Part 1 of 
the figure shows a permutation tt with a long cycle on the torus which almost meets itself in the x direction. 
In part 2, after a Metropolis step sending tt to tt' , one cycle winds by +1 and the other by —1. Metropolis 
steps create winding cycles only in opposite-direction pairs; the total W x (tt) is still zero. Part 3 of the figure 
shows that if we reverse one cycle (which is a zero-energy move), W x (tt) is now 2. In general (with full details 
in |Kerl] ) . winding numbers of even parity can be generated. 

Our current best algorithm (swap-and-reverse) has two types of sweeps: (1) For each lattice site, do a 
Metropolis step as above. (2) For each cycle in the permutation, reverse the direction of the cycle with 
probability 1/2. This permits winding numbers of even parity in each of the three axes. 

We have experimented with various methods to obtain winding numbers of all parities. The creation or 
destruction of a winding cycle is a non-local update; one is reminded of the Swendsen-Wang algorithm for 
the Ising model. However, our attempt at non-local updates has an unreasonably low acceptance rate, 
namely, on the order of e~ L where L is the box length. 
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Figure 7: Conservation of winding number in the swap-only algorithm. 



We have also created a worm algorithm, inspired by approaches to this same winding-number problem in 
path-integral Monte Carlo methods [BPS06, PST98 . That is, a permutation loop is selected at random and 
then cut open at a randomly selected point. The resulting worm is allowed to move around A via Metropolis 
moves; eventually it closes again. This worm algorithm has an elegant theory and correctness proof [Kerlj ; 
yet, it has an unacceptably long stopping time for loop closure, and none of our attempts to remedy the 
stopping-time problem have satisfied detailed balance. 

At present, we content ourselves with the swap-and-reverse algorithm; it is used to generate all the results 
discussed in section^] The order parameters fs and fw depend on winding phenomena, but the other three, 
l/£, //, and / max , do not; furthermore, results obtained in section 0] using each of the five order parameters 
are, for the most part, compatible. Yet, as we will see, fs and fw do not permit successful finite-size scaling. 

3.5 Finite-size scaling 

Finite-size scaling takes the form of a hypothesis, or rather a set of hypotheses, which is tested against the 
data. See also [CG GPj for a nice survey. 

We have an infinite- volume random variable S(T), e.g. any of the order parameters defined in section [2~3l 
The finite-volume quantity is S L (T). Define t = (T — T c )/T c . Examine, say, 0.99 < t < 1.01. The first 
hypothesis is that the correlation length £(T) follows a power law 

t;(T)~\t\-», T^T C 

For the infinite-volume quantity, we expect a power-law behavior 

S(T)~f, (-*)', or \t\". 

(The domain of validity is t < or t > depending on whether the order parameter S is left-sided or 
right-sided, respectively.) One moreover hypothesizes that for T near T c , Sl(T) and S(T) are related by a 
universal function Qs which depends on T only through the ratio L/£: 

S L (T) = L~ p/ "Q s (L 1/v t) ~ L-o/»Q s ((L/0 1/u )- (3-13) 



4 Results 

Here we complete the steps sketched in section l3~5l The flow of data and uncertainties are as follows: 

• Markov chain Monte Carlo simulations, with error bars determined using the method of integrated 
autocorrelation time [Berg], yield Sl(T,cx) data points. There are five order parameters 5, six values 
of L (30, 40, 50, 60, 70, 80), nine values of a, and a few dozen values of T for each a. 
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/max' = 




/ msD[ , a = 0.004 




Figure 8: Order parameters /m and l/£ for L = 40,60,80 and a = and 0.004. The remaining order 
parameters fg, fw, and // behave similarly to Jm but with not all with the same critical exponents. 



• CPU time per L,T,a experiment, with 10 5 Metropolis sweeps, is approximately 1.3 hours for L = 40 
and 20 hours for L = 80. For the work described in this paper and in Kerl , a total of 5.4 CPU years 
was used. 

• For each S, L, and a, we use Sx(T, a) values for all available values of T and a to estimate pg(L). 
(Critical exponents are assumed to be independent of a for small a, or with weak enough dependence 
on a that that dependence is lost in the noise.) Error bars may be propagated from the MCMC 
simulations, or computed from regression uncertainties. 

• Extrapolating ps{L) in L — > oo results in the five estimated critical exponents pg. Uncertainties are 
computed from the regression analysis. 

• Once the critical exponents are estimated, we obtain T Cj s(a) for each of the five order parameters S 
and for each a. Uncertainties are computed by visual inspection of the crossing plots discussed in 
section 14.31 

• Once the critical exponents and T c are known, one should be able to obtain plots of the universal 
function Qg which is, up to sampling variability, independent of L, T, and a. This verifies the finite- 
size-scaling hypothesis. 

• The shift in reduced critical temperature is as in equation (|2.10p . Error bars are computed from 
regression uncertainties. 



4.1 Determination of L-dependent critical exponents 

For each of order parameter S, interaction parameter a, and box length L, we examine all S(L, T, a) data for 
which S > s, with e taken from the plots to ensure that we examine the portions of the curves corresponding 
to non-zero order parameter in the infinite limit (see figure [HJ . For l/£, this means T > T c ; for the other four 
order parameters, this means T < T c , From plots such as those in figure we choose e to be 0.1 for l/£, 
0.01 for Jm, 0.01 for //, 0.05 for fg, and 0.01 for fw For each S, a, and L, we then apply linear regression 
to S'(L,T) 1 / ps for varying pg. We find ps{L) which optimizes the correlation coefficient | Young | of the 
linear regression. Results are shown in figure Given ps{L) along with its corresponding linear-regression 
parameters m and b, we may plot a power-law fit to the simulational data. One such comparison plot is 
shown in figure ITOl 
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Figure 9: On the left: determination of critical exponent ps(L, a) for order parameter as the value which 
minimizes linear-regression error for SL{T,ay/P. Visually, one sees ps{L = 80, a = 0.0) « 0.59. On the 
right: estimated critical exponents for L = 30, 40, 50, 60, 70, 80. 



a 


Mean 


Std.err. 


Count 


0.000 


0.6242981 


0.0000897 


78 


0.0001 


0.6243312 


0.0001079 


78 


0.0002 


0.6245691 


0.0000921 


72 


0.0005 


0.6245402 


0.0001062 


66 


0.0008 


0.6244347 


0.0000856 


72 


0.001 


0.6244779 


0.0001020 


60 


0.002 


0.6246345 


0.0001154 


42 


0.003 


0.6245906 


0.0001559 


48 


0.004 


0.6245966 


0.0001964 


42 



Table 1: /m/ fi as a function of a. An upward trend is visible, but it is not pronounced. 




V 


0.5559 


± 0.0037 


Ps 


0.6201 


± 0.0065 


Pw 


0.7750 


± 0.0073 


Pi 


0.7451 


± 0.0052 


Pm 


0.7486 


± 0.0059 



Table 2: Extrapolated estimates of the infinite- volume critical exponents, found from the vertical intercept 
of figure [5] 




6.90 



Figure 11: The crossing method to estimate T c (a) for order parameter //, with p and v as above: T c (a) 
corresponds to the horizontal coordinate of the intersection point of the plots. The upper-right-hand plot is 
a close-up of the upper-left- hand plot. Order parameters fs and fw, which depend on winding phenomena, 
do not exhibit clear crossing behavior. 



4.2 Extrapolation of critical exponents for the infinite-volume limit 

Next, for each S, given estimates ps(L) for increasing values of L, we plot ps(L) versus 1/L. The vertical 
intercept of this plot estimates the infinite- volume exponent ps{a). (See figure 0) Results are shown in 
tabled 



4.3 Determination of critical temperature 

Given the above estimators of the critical exponents, the crossing method [CGGPj estimates T c (a). Namely, 
we plot LP/*S L (T) as a function of T. At T = T c we have t = and Lp/»S l {T) = Q s (0), regardless of L 
(equation (|3.13jl ). Thus, these curves will cross (approximately, due to sampling variability) at T = T c . If 
they do not, the finite-size-scaling hypothesis is not verified. (Note in particular that for order parameter 
l/£ whose critical exponent is u, we apply the crossing method to LSl(T) as a function of T: thus, the T c (a) 
estimate using l/£ is independent of v.) See for example figure [TTJ (We acknowledge that larger values of 
L, would improve the visual effect. Results presented here are those obtained within the timeframe of the 
author's doctoral dissertation work.) Results are in figure H"3l 

Using order parameters fs and fw, which depend on winding phenomena, one does not see clear crossing 
behavior. We suggest that either this is related to the even-winding-number issue discussed in section \3A\ 
or fs and fw are not good order parameters for this model. We suspect the former; in every manner except 
this crossing issue, fs and fw behave as expected. (In the absence of clear crossing behavior for fs and fw, 
for the sake of discussion we nonetheless provide best visual estimates for T c (a) for fs and fw- These will 
not be used for further analysis toward our final result.) 
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Figure 13: Critical temperature as function of a. 



4.4 Verification of finite-size-scaling hypothesis 

Now that we have estimated ps, v, and T c (a) for each of the five order parameters 5, we may plot 
a) as a function of L x l v t. This is a plot of the scaling function Qg- If the hypothesis is correct, 
the curves for all L should coincide, or collapse, to within sampling error — which they do (e.g. figure IT2"j) . 



4.5 Determination of the shift in critical temperature 

As discussed in section |2~4"[ we are seeking a linear relationship between AT c (a) and a, with constant c. This 
can be visualized in figure IM) which is obtained from the T Cj s(a) data of figure [13] using equation (|2.10[) . 
We start with all the (a, AT c (a)) data points from section l4~3l We omit values obtained using fs and fw, 
due to the aforemention lack of crossing behavior. We also omit values obtained using a = 0.004, since the 
critical-temperature plots of figure [131 suggests that this starts to exceed the domain of linear approximation. 
We perform a linear regression with error bars [Young| on the (a, AT c (a)) data points. We use a slope-only 
fit, rather than a slope-intercept fit, since AT c {a) has zero intercept by its definition. We find 

c = 0.618 ± 0.086 (2 a error bar). 

Within experimental uncertainty, this result, for points on the lattice with Ewens cycle- weights, matches the c 
value of equation (|2.1ip for point positions varying on the continuum with decaying-cycle- weight interactions. 
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Figure 14: Shift in critical temperature, and linear fit, as function of a. Recall from equation (12.101) that 
AT c (a) = jTq^ • Order parameters fg and fw were omitted from the fit, due to lack of crossing 
behavior; a = 0.004 was omitted due to onset of curvature of T c (a). The heavy solid line shows a linear 
fit with empirically determined constant of proportionality; the lighter solid line is the comparison value of 
Betz and Ueltschi (slope 2/3) for decaying cycle weights and continuum point positions. 



4.6 Constancy of the macroscopic-cycle quotient 

As discussed in section |2~4"1 we hypothesize that the macroscopic-cycle quotient fu/fi in the infinite- volume 
limit is dependent on a but is constant in T where it is defined, i.e. for T < T c since fi = for T > T c . 
This may be visualized by comparing figures such as El one sees that /m and fi appear to have the same 
critical exponent. Alternatively, one may plot the ratio fu/fi f figure [15]) . In the infinite-volume limit, fx is 
zero for T > T c and so we are interested only in the values of the quotient for T < T c . In that region, the 
quotient does indeed appear to be constant in T. 

We test this constancy hypothesis as follows. The respective critical exponents are pu and pj. The estimators 
are pm and pi, computed by averaging over several different values of L and a as described in section l4~2l 
Treating these estimators as normally distributed (as justified by the raw data), we obtain the standard 
deviations of the pm,i(L, a) samples, along with the standard deviations of the means pm,i- 

0.7482 Pl = 0.7445 

0.0428 si = 0.0374 

50 m = 50 

0.006059 s//Vn7 = 0.005295. 

The difference pm — pi is also normally distributed about the true mean pm — pi, but pm and pi are not 
independent since they are sample means of random variables computed from the same Markov chain Monte 
Carlo sequence of permutations. Thus 

V&r(p M - pi) = Var(p M ) + Var( j o / ) - 2Cov(p M ,Pi)- 

Computing the sample covariance of the pi\i(L,a) and pi(L,a) data series, we obtain the covariance and 
resulting standard error Sd of the difference 

Cov(p M ,Pi) = 0.0004 s d /VH = 0.0070. 
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Figure 15: Macroscopic-cycle quotient /m//j for a = 0,0.002. 



Normalizing, we find 



pit- Pi = 0-0037 



Pm - Pi 



0.0037 



Sdj^n 0.0070 



= 0.5293. 



We hypothesize Pm ~ Pi = 0; the estimated value Pm ~ Pi lies comfortably within a standard deviation of 
this. We note, moreover, that the value of /m//j, while constant in T, trends upward with a (see table [T] 
and figure [T6|) . This merits further investigation. 



4.7 Conclusions 



(1) For annealed point positions, equation {HT3) gives T c (0) w 6.625. Our result T c (0) = 6.873 ± 0.006 
(2tr error bar) unambiguously shows that the lattice structure modifies the critical temperature, even in the 



, +6.241e-l 



/max/// as function of a 




0.004 



Figure 16: /m//j as a function of a. 
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non-interacting (a = 0) case. 



(2) As detailed in section |4~5I we find that the reduced shift in critical temperature as a function of interaction 
parameter a is 

. . , T c (a)-T c (0) 
V ' T c (0) 

with 

c = 0.618 ±0.086 (2ct error bar). 

This is compatible (section [2 .4[) with the related result of [BU08] . Even though the lattice structure changes 
the critical temperature (conclusion 1), the shift in critical temperature is unaffected. 

(3) As described in section [2T4l Shepp and Lloyd |SLj find that E[£ max ]/A ~ 0.6243 for uniform-random 
(non-spatial) permutations. For spatial permutations, we define a macroscopic-cycle quotient E[^ roax ]/JV/i 
which is the ratio of mean maximum cycle length as a fraction of the number of sites in long cycles. Our 
result (table [1]) is compatible with that of Shepp and Lloyd for the non-interacting case, with an increase 
which appears to be linear as a function of interaction parameter a. Our result is also compatible with 

GRU , which (among other conclusions) recovered the Shepp and Lloyd result for the a = case. 



5 Future work 



Now that the a-dependence of the macroscopic-cycle quotient's constant upon a has been found empirically, 
one would next like to explain that dependence analytically. 

Ideally, one would have an algorithm to permit odd winding numbers, as discussed in section [3. 41 

Sampling from the true Bose-gas distribution using the random-cycle model requires three changes. First, one 
needs to conduct simulations using the Bose-gas interaction (equation (I2.ip ) rather than the cycle- weight 
interaction (equation (12.21) ). The interaction term V is a CPU- intensive Brownian-bridge computation 
BU07 ; unpublished work of Ueltschi and Betz shows that it may be approximated in the weak- interaction 
case by a simpler Riemann integral. Precomputed tables and interpolation may make use of this integral 
feasible. Second, point positions must be allowed to vary on the continuum. This entails a second type of 
Metropolis step, in addition to that shown in section I3TT1 Third, since points are no longer held fixed on the 
lattice, it is no longer trivial to find nearest neighbors. Software efficiency requires a hierarchical partitioning 
of A. The second and third points simply require a software effort. Implementing them will be worthwhile 
only if the interaction terms can be simplified to the point that they are computationally feasible, which is 
a mathematical effort. 
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